Developmental and Environmental Influences on Thalamocortical Connectivity are Patterned Along the S-A Axis in a Clinically Enriched Sample
Read in Data for Developmental and Environmental Characterization
Glasser parcel names/tracts
#Glasser regions with corresponding labels and tract names
glasser.regions <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360_regionlist.csv") #parcel name, label name
glasser.regions$tract <- paste0("thalamus_", glasser.regions$orig_parcelname) %>% gsub("_ROI", "_autotrack", .) %>% gsub("-", "_", .) S-A axis
S.A.axis.cifti <- read_cifti("/cbica/projects/thalamocortical_development//Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo
S.A.axis <- data.frame(SA.axis = rank(S.A.axis.cifti$data), orig_parcelname = names(S.A.axis.cifti$Parcel))
S.A.axis <- merge(S.A.axis, glasser.regions, by = c("orig_parcelname"), sort = F)
S.A.axis$SA.axis.bin <- as.numeric(cut2(S.A.axis$SA.axis, g = 5)) #divide the S-A axis into 5 ranked bins, 72 regions per binNeuroSynth cognitive atlas term maps
#Neurosynth z-score maps generated by nimare. Nimare uses multilevel kernel density analysis- Chi-square analysis + FDR-correction to use the same procedure as Neurosynth
neurosynth.terms <- read.csv("/cbica/projects/thalamocortical_development/Maps/neurosynth/atl-glasser360_neurosynth.csv") %>% select(-regionName) %>% select(-timing) #neurosynth term meta-analytic scores for 123 terms present in both NeuroSynth and the Cognitive Atlas, ordered in lh --> rh glasser order
neurosynth.termlist <- names(neurosynth.terms)
neurosynth.terms$label[1:180] <- glasser.regions$label[181:360] #lh first
neurosynth.terms$label[181:360] <- glasser.regions$label[1:180] #rhSpatial permutation (spin) test parcel rotation matrix
perm.id.full <- readRDS("/cbica/projects/thalamocortical_development/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds") #10,000 spatial null spins
spin.parcels <- glasser.regions #order of complete set of glasser parcels for spinningDataset-specific tract list
tractlist.hbn <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HBN/HBN_tractlist_primary.csv") #generated by tract_measures/tractlists/thalamocortical_tractlists.RDevelopmental statistics
setwd("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/") #gam output dir
#list all FA development files in development results directory for accessing
files <- list.files(getwd(), pattern='FA', ignore.case=T, full.names = F)
files <- files[files %>% grep("hcpd", ., invert = T)] #only files for HBN and (to assess replication) PNC
filenames <- c()
#generate variable names to assign files data to
for(name in files){
Rname <- gsub('^.{12}|.{4}$', '', name) #remove ".csv" from end of filename and "development_" from start of filename
filenames <- append(filenames, Rname) #save filenames into a character vector
}
#read in files and assign to variables
for(i in 1:length(filenames)){
Rfilename <- sprintf("%s",filenames[i]) #save filename as an individual string
if(grepl("gameffects", Rfilename)) {
x <- read.csv(files[i], header = TRUE)
x <- merge(x, glasser.regions, by = "tract")
assign(Rfilename, x)
}
if(!grepl("temporal", Rfilename) & !grepl("gameffects", Rfilename)) {
x <- read.csv(files[i], header = TRUE)
x <- merge(x, S.A.axis, by = "tract")
assign(Rfilename, x)
}
}
rm(x)Environment statistics
setwd("/cbica/projects/thalamocortical_development/thalamocortical_results/environment_results/") #gam output dir, from gam_functions/fit_envGams.R
#list all .csv files in the environment results directory
files <- list.files(getwd(), pattern = '.csv', ignore.case = T, full.names = F)
files <- files[files %>% grep("hbn", ., invert = F)] #only files for HBN and (to assess replication) PNC
for(i in 1:length(files)){
filename <- files[i]
Rname <- gsub('.{4}$', '', filename)
x <- read.csv(filename, header = TRUE)
assign(Rname, x)
}
rm(x)A Spectrum of Thalamocortical Developmental Trajectories
Cortex-wide thalamic connectivity developmental profiles
Number of significant developmental effects
#Function to calculate the number and percent of thalamocortical connections showing a significant developmental change in a given measure
sig.effects <- function(measure, atlas, dataset){
ageeffects.df <- get(sprintf("gameffects_%s_%s_%s", measure, atlas, dataset))
ageeffects.df$significant <- p.adjust(ageeffects.df$Anova.smooth.pvalue, method = c("fdr")) < 0.05 #fdr-corrected significant connections
sigeffect.totaln <- sum(ageeffects.df$significant) #total number of significant connections
sigeffect.percent <- round(sigeffect.totaln/length(ageeffects.df$significant)*100, 2) #percent of significant connections
print(sprintf("In %s, %s thalamocortical connections (%s percent) show significant age-related changes in %s", dataset, sigeffect.totaln, sigeffect.percent, measure))
}HBN
## [1] "In hbn, 166 thalamocortical connections (73.78 percent) show significant age-related changes in FA"
Thalamocortical connection GAM smooth functions
smoothcentered_FA_glasser_hbn.plot <- merge((smoothcentered_FA_glasser_hbn %>% filter(grepl("thalamus_R", tract))), (gameffects_FA_glasser_hbn %>% select(tract, GAM.smooth.partialR2)), by = "tract") #zero-centered developmental smooth functions for RH connections + tract-specific partial R2 for plotting
ggplot(smoothcentered_FA_glasser_hbn.plot, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) +
geom_line(linewidth = .3, alpha = .8) +
scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.017, 0.075), oob = squish) +
theme_minimal() +
labs(x="\nAge", y="Connection FA\n") +
scale_x_continuous(breaks=c(6, 8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
theme(
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))Maturational patterns are highly reproducible
Correlation of thalamocortical connection development effect sizes between datasets
Pearson’s R
gameffects.merged <- merge(gameffects_FA_glasser_pnc, gameffects_FA_glasser_hbn, by = c("tract", "label", "orig_parcelname"), suffixes = c("_pnc", "_hbn"))
cor.test(gameffects.merged$GAM.smooth.partialR2_pnc, gameffects.merged$GAM.smooth.partialR2_hbn)##
## Pearson's product-moment correlation
##
## data: gameffects.merged$GAM.smooth.partialR2_pnc and gameffects.merged$GAM.smooth.partialR2_hbn
## t = 15.616, df = 218, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6575492 0.7835828
## sample estimates:
## cor
## 0.7266226
Spatial permutation (spin) based p-value
spin.data <- left_join(spin.parcels, gameffects.merged, by = c("label", "orig_parcelname", "tract")) #full set of parcel data in rh --> lh order for spinning. spin test null correlations use complete obs only. Each null correlation correspondence statistic is thus computed on a slightly reduced set of data, akin to a jackknife procedure
perm.sphere.p(x = spin.data$GAM.smooth.partialR2_pnc, y = spin.data$GAM.smooth.partialR2_hbn, perm.id = perm.id.full, corr.type = "pearson") ## [1] 5e-05
Correlation plot
ggplot(gameffects.merged, aes(x = GAM.smooth.partialR2_pnc, y = GAM.smooth.partialR2_hbn)) +
geom_point(color = c("#FCAB6A"), size = 0.5) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
theme_classic() +
scale_x_continuous(limits = c(-0.01, 0.17)) +
scale_y_continuous(limits = c(-0.051, 0.15)) +
labs(x="\nPNC", y="HBN\n") +
theme(
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/FA_Ageeffect_correlation_PNCHBN.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.6)Correlation of thalamocortical connection age of maturation between datasets
Pearson’s R
cor.test(gameffects.merged$smooth.increase.offset_pnc, gameffects.merged$smooth.increase.offset_hbn)##
## Pearson's product-moment correlation
##
## data: gameffects.merged$smooth.increase.offset_pnc and gameffects.merged$smooth.increase.offset_hbn
## t = 6.029, df = 151, p-value = 1.214e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3029697 0.5600095
## sample estimates:
## cor
## 0.4404722
Spatial permutation (spin) based p-value
perm.sphere.p(x = spin.data$smooth.increase.offset_pnc, y = spin.data$smooth.increase.offset_hbn, perm.id = perm.id.full, corr.type = "pearson") ## [1] 0.0015
Early and Late Maturing Thalamocortical Connections
NeuroSynth decoding of age of maturation maps
Compute cognitive term-development map spatial correlations
#Function to calculate the correlation between a neurosynth term z-score map and a thalamocortical development map
developmentmap.neurosynth.decoding <- function(df, dev.metric, term){
df.allregions <- left_join(spin.parcels, df, by = c("label", "tract", "orig_parcelname")) #full set of parcel data in rh --> lh order
term.cor <- cor(df.allregions[,term], df.allregions[,dev.metric], use = "complete.obs", method = c("spearman")) #correlation between neurosynth term map and development metric
term.results <- data.frame("term" = term, "term.correlation" = term.cor)
return(term.results)
}Developmental timing decoding by Neurosynth cognitive terms
HBN
gameffects_neurosynth_hbn <- merge(gameffects_FA_glasser_hbn, neurosynth.terms, by = "label")
devpattern.neurosynth.hbn <- ldply(neurosynth.termlist, function(t){
developmentmap.neurosynth.decoding(df = gameffects_neurosynth_hbn, dev.metric = "smooth.increase.offset", term = t)}) %>% #neurosynth-neurodevelopment correlations for all terms
arrange(desc(term.correlation))devpattern.neurosynth.hbn <- rbind(slice_max(devpattern.neurosynth.hbn, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.hbn, order_by = term.correlation, n = 10)) #10 most positively correlated and negatively correlated neurosynth terms with the age of maturation map
devpattern.neurosynth.hbn## term term.correlation
## 1 cognitive_control 0.6111711
## 2 recall 0.5162633
## 3 thought 0.5000909
## 4 decision 0.4826017
## 5 expectancy 0.4703547
## 6 memory 0.4511850
## 7 retrieval 0.4442939
## 8 memory_retrieval 0.4397963
## 9 context 0.4276347
## 10 judgment 0.4230770
## 11 multisensory -0.5788485
## 12 visual_perception -0.5575273
## 13 object_recognition -0.5412269
## 14 gaze -0.5204399
## 15 localization -0.4933611
## 16 detection -0.4488658
## 17 perception -0.4268284
## 18 facial_expression -0.3828248
## 19 movement -0.3641246
## 20 visual_attention -0.3538157
ggplot(devpattern.neurosynth.hbn, aes(x = term.correlation, y = term, color = term.correlation)) +
geom_segment(aes(y = reorder(term, term.correlation), yend = term, x = 0, xend =
term.correlation), color = "#989898", linewidth = 0.2) +
geom_point(size = 1.8, alpha = 0.85) +
theme_light() +
scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.29, 0.35), oob=squish) +
labs(x = "\nCognitive Term-Thalamocortical Development Correlation", y = "NeuroSynth Term\n") +
scale_x_continuous(limits = c(-0.6,0.62), breaks = c(-0.6, -0.3, 0, 0.3, 0.6)) +
theme(
legend.position = "none",
panel.grid.major.y = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.y = element_blank(),
panel.grid.minor.x = element_blank(),
axis.ticks.y = element_blank(),
panel.border = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.ticks = element_line(linewidth = .2))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/Neurosynth_Maturationdecoding_HBN.pdf", device = "pdf", dpi = 500, width = 1.9, height = 3.1)Overlapping terms across datasets
gameffects_neurosynth_pnc <- merge(gameffects_FA_glasser_pnc, neurosynth.terms, by = "label")
devpattern.neurosynth.pnc <- ldply(neurosynth.termlist, function(t){
developmentmap.neurosynth.decoding(df = gameffects_neurosynth_pnc, dev.metric = "smooth.increase.offset", term = t)}) %>% #neurosynth-neurodevelopment correlations for all terms
arrange(desc(term.correlation))
devpattern.neurosynth.pnc <- rbind(slice_max(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10)) #10 most positively correlated and negatively correlated neurosynth terms with the age of maturation mapmerge(devpattern.neurosynth.pnc, devpattern.neurosynth.hbn, by = "term") %>% arrange(desc(term.correlation.x)) %>% select(term)## term
## 1 expectancy
## 2 cognitive_control
## 3 memory_retrieval
## 4 thought
## 5 memory
## 6 recall
## 7 facial_expression
## 8 perception
## 9 gaze
## 10 visual_perception
## 11 object_recognition
Thalamocortical Connectivity Maturation is Organized by the Sensorimotor-Association Axis
gameffects_FA_glasser_hbn <- merge(gameffects_FA_glasser_hbn, (S.A.axis %>% select(tract, SA.axis)), by = "tract")Thalamocortical connections mature at progressively older ages across the S-A axis
HBN
derivatives_FA_glasser_hbn$significant.derivative[derivatives_FA_glasser_hbn$significant.derivative <= 0] <- NA #replace non-significantly increasing derivatives with NA. (Note here 0 was == not significant)
ggplot() +
geom_line(data = derivatives_FA_glasser_hbn, aes(x = age, y = SA.axis, group = tract, alpha = significant.derivative, color = SA.axis, linewidth = significant.derivative)) +
scale_alpha_continuous(range = c(0, 1), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
scale_linewidth_continuous(range = c(0, 4.5), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
scale_color_gradient2(low = "#ffc933", mid = "#e9dcf2", high = "#6f1282", guide = "colorbar", midpoint = 180) +
scale_x_continuous(breaks = c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .25)) +
scale_y_continuous(expand = c(0.03, 0)) +
theme_classic() +
ylab("Sensorimotor-Association Axis Rank\n") +
xlab("\nAge") +
theme(
legend.position = "none",
axis.line = element_line(linewidth = .2),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.ticks = element_line(linewidth = .2))Spearman’s rho
cor.test(gameffects_FA_glasser_hbn$SA.axis, gameffects_FA_glasser_hbn$smooth.increase.offset, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: gameffects_FA_glasser_hbn$SA.axis and gameffects_FA_glasser_hbn$smooth.increase.offset
## S = 196170, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.6899525
Spatial permutation (spin) based p-value
spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hbn, by = c("label", "orig_parcelname", "tract"))
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") ## [1] 0.00205
Correlation plot
ggplot(gameffects_FA_glasser_hbn, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
geom_point(size = 0.5) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
theme_classic() +
scale_y_continuous(limits = c(11, 20), breaks = c(12, 14, 16, 18, 20)) +
labs(x="\nS-A axis", y="Age of maturation (years)\n") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) Thalamocortical Connectivity Maturation Mirrors Timescales of Cortical Plasticity
Plasticity marker development data
#Fluctuation amplitude age of decrease onset data from Sydnor et al., 2023, Nat Neurosci https://www.nature.com/articles/s41593-023-01282-y
boldamplitude.agedecrease.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/spatiotemp_dev_plasticity/cortical_maps/AgeofDeclineOnset_FirstNegDeriv.pscalar.nii")
boldamplitude.dev <- data.frame(orig_parcelname = names(boldamplitude.agedecrease.cifti$Parcel), amplitude.agedecline = boldamplitude.agedecrease.cifti$data)#T1/T2 ratio rate of developmental change from Baum et al., 2022, J Neuro https://www.jneurosci.org/content/42/29/5681
myelin.maxdev.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_median_posterior_age_of_max_slope_myelination.pscalar.nii") #age of maximal T1/T2 ratio increase
myelin.roc.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_mean_posterior_annualized_roc_myelin.pscalar.nii") #annualized rate of change
myelin.dev <- data.frame(orig_parcelname = names(myelin.maxdev.cifti$Parcel), myelin.agemaxdev = myelin.maxdev.cifti$data, myelin.annualroc = myelin.roc.cifti$data)Thalamocortical maturational patterns align to the spatiotemporal unfolding of plasticity readouts
HBN
plasticity.dev.hbn <- left_join(plasticity.dev, gameffects_FA_glasser_hbn, by = c("label", "orig_parcelname", "tract"))Thalamocortical connectivity maturation - E:I ratio magnitude of developmental decline
Spearman’s rho
cor.test(plasticity.dev.hbn$smooth.increase.offset, plasticity.dev.hbn$EI.ageslope, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: plasticity.dev.hbn$smooth.increase.offset and plasticity.dev.hbn$EI.ageslope
## S = 272148, p-value = 8.219e-15
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.5698692
Spatial permutation (spin) based p-value
EI.thalamus.hbn <- perm.sphere.p(x = plasticity.dev.hbn$smooth.increase.offset, y = plasticity.dev.hbn$EI.ageslope, perm.id = perm.id.full, corr.type = "spearman") Correlation plot
ggplot(plasticity.dev.hbn, aes(x = EI.ageslope, y = smooth.increase.offset)) +
geom_point(size = 0.5, color = "#fcd16d") +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_y_continuous(limits = c(12, 18), breaks = c(12, 14, 16, 18)) +
theme_classic() +
labs(x="E/I ratio decline (age slope)\n", y="\nThalamocortical age of maturation (years)") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) Thalamocortical connectivity maturation - myelin annualized rate of growth
Spearman’s rho
cor.test(plasticity.dev.hbn$smooth.increase.offset, plasticity.dev.hbn$myelin.annualroc, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: plasticity.dev.hbn$smooth.increase.offset and plasticity.dev.hbn$myelin.annualroc
## S = 998646, p-value = 2.643e-15
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.5783632
Spatial permutation (spin) based p-value
myelin.thalamus.hbn <- perm.sphere.p(x = plasticity.dev.hbn$smooth.increase.offset, y = plasticity.dev.hbn$myelin.annualroc, perm.id = perm.id.full, corr.type = "spearman") Correlation plot
ggplot(plasticity.dev.hbn, aes(x = myelin.annualroc, y = smooth.increase.offset, color = myelin.annualroc)) +
geom_point(size = 0.5, color = "#F1A45F") +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_x_continuous(limits = c(0.005, 0.025)) +
scale_y_continuous(limits = c(12, 18), breaks = c(12, 14, 16, 18)) +
theme_classic() +
labs(x="T1/T2 ratio annualized growth rate\n", y="\nThalamocortical age of maturation (years)", ) +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) Thalamocortical connectivity maturation - fluctuation amplitude age of decrease onset
Spearman’s rho
cor.test(plasticity.dev.hbn$smooth.increase.offset, plasticity.dev.hbn$amplitude.agedecline, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: plasticity.dev.hbn$smooth.increase.offset and plasticity.dev.hbn$amplitude.agedecline
## S = 147055, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.6852288
Spatial permutation (spin) based p-value
amplitude.thalamus.hbn <- perm.sphere.p(x = plasticity.dev.hbn$smooth.increase.offset, y = plasticity.dev.hbn$amplitude.agedecline, perm.id = perm.id.full, corr.type = "spearman")Correlation plot
ggplot(plasticity.dev.hbn, aes(x = amplitude.agedecline, y = smooth.increase.offset, color = amplitude.agedecline)) +
geom_point(size = 0.5, color = "#7B2C80") +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_x_continuous(limits = c(8.2, 21.2), breaks = c(8, 10, 12, 14, 16, 18, 20)) +
scale_y_continuous(limits = c(12, 18), breaks = c(12, 14, 16, 18)) +
theme_classic() +
labs(x="Amplitude decrease onset (years)\n", y="\nThalamocortical age of maturation (years)") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) FDR-corrected p-values across correlations
## [1] 0.00280 0.01495 0.00450
## [1] 0.00675 0.01495 0.00675
Lolly plot of correlation strength
plasticity.results <- data.frame(map = factor(), corr = double())
plasticity.results <- plasticity.results %>% add_row(map = "E/I Ratio", corr = 0.5698692)
plasticity.results <- plasticity.results %>% add_row(map = "T1/T2 Ratio", corr = -0.5783632)
plasticity.results <- plasticity.results %>% add_row(map = "BOLD Amplitude", corr = 0.6852288)
plasticity.results$map <- factor(plasticity.results$map, ordered = TRUE, levels = c( "E/I Ratio", "T1/T2 Ratio", "BOLD Amplitude"))
ggplot(plasticity.results, aes(x = map, y = corr)) +
geom_segment(aes(x = map, xend = map, yend =
corr, y = 0), linewidth = 0.2, color = "#989898") +
geom_point(size = 2, alpha = 0.85, color = "#FCAB6A", fill = "#FCAB6A", shape = 22) +
labs(x="Cortical Development Map") +
labs(y="\nThalamocortical Correlation (r)") +
theme_classic() +
scale_y_continuous(breaks = c(-0.6, -0.3, 0, 0.3, 0.6), limits = c(-0.62, 0.7)) +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.text.x = element_text(vjust = 0.7, angle = 20),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))Associations between Neighborhood Socioeconomic Conditions and Thalamocortical Connectivity
Number of significant environmental effects
Significant neighborhood-level environment effects (neighborhood SES)
generalSES_maineffects_FA_glasser_hbn$significant <- p.adjust(generalSES_maineffects_FA_glasser_hbn$GAM.cov.pvalue, method = c("fdr")) < 0.05
sum(generalSES_maineffects_FA_glasser_hbn$significant)## [1] 120
sum(generalSES_maineffects_FA_glasser_hbn$significant)/length(generalSES_maineffects_FA_glasser_hbn$significant)## [1] 0.5333333
Extract significant connections
Percent of positive significant neighborhood-level environment effects
## [1] 100
Neighborhood Environment Effects Vary across the Sensorimotor-Association Axis
S-A axis correlation
Spearman’s rho
spin.data <- left_join(spin.parcels, sigeffects.df, by = "tract") #merged data for spin tests, with all 360 parcels in the expected right hemisphere --> left hemisphere order (matching glasser.coords)
spin.data <- merge(spin.data, S.A.axis, by = c("label", "orig_parcelname", "tract"), sort =F)
cor.test(spin.data$GAM.cov.tvalue, spin.data$SA.axis, method = c("spearman"))##
## Spearman's rank correlation rho
##
## data: spin.data$GAM.cov.tvalue and spin.data$SA.axis
## S = 199426, p-value = 0.0006705
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.3075005
Spatial permutation (spin) based p-value
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$GAM.cov.tvalue, perm.id = perm.id.full, corr.type = "spearman") ## [1] 0.0349
Correlation plot
ggplot() +
geom_point(data = spin.data, aes(x = SA.axis, y = GAM.cov.tvalue, color = SA.axis), size = 0.5) +
geom_smooth(data = spin.data, aes(x = SA.axis, y = GAM.cov.tvalue), method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
theme_classic() +
labs(x="\nS-A axis", y="Environment effect (t)") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) Effect magnitude enrichment testing (S-A axis quintiles)
Mean environment effect t-value by S-A axis bin (empirical)
#Calculate empirical enrichment statistic for each S-A axis bin
envSES_effects_SAaxisbins <- spin.data %>%
group_by(SA.axis.bin) %>%
do(mean.tvalue.significant = mean(.$GAM.cov.tvalue, na.rm = TRUE)) %>%
unnest(cols = c("mean.tvalue.significant"))## # A tibble: 5 × 2
## SA.axis.bin mean.tvalue.significant
## <dbl> <dbl>
## 1 1 2.84
## 2 2 2.88
## 3 3 2.92
## 4 4 2.97
## 5 5 3.22
Null distribution and enrichment testing by S-A axis quintile
#Calculate null distribution of t-values by bin based on spherical spatial permutations (spins)
## inputs to spatially permute (i.e., spin) with the pre-computed spatial permutation matrix
x <- spin.data$SA.axis.bin #cortical map 1 to spatially permute: S-A axis bins
y <- spin.data$GAM.cov.tvalue #cortical map 2 to spatially permute: significant envSES-thalamocortical FA t-values
perm.id <- perm.id.full
## spin the empirical data
nroi = dim(perm.id)[1] #number of regions
nperm = dim(perm.id)[2] #number of permutations
x.perm = y.perm = array(NA,dim=c(nroi,nperm)) #initialize
for (r in 1:nperm) {
for (i in 1:nroi) {
x.perm[i,r] = x[perm.id[i,r]] #spinning x, spatially permuted bins
y.perm[i,r] = y[perm.id[i,r]] #spinning y, spatially permuted t-values
}
}
## compute mean t-value for each S-A axis bin using spatially permuted x inputs
### set up spun (x) and empirical (y) df
x.spatialperm <- as.data.frame(x.perm) %>% set_names(sprintf("perm%s", seq(from = 1, to = ncol(x.perm))))
x.spatialperm.empiricaly <- cbind(y, x.spatialperm) %>% as.data.frame()
colnames(x.spatialperm.empiricaly)[1] <- c("empirical.y")
### compute mean t-values for S-A axis bins for all spins
envSES_effects_SAaxisbins_permutedx <- sapply(x.spatialperm.empiricaly[2:ncol(x.spatialperm.empiricaly)], #for all permuted bins
function(p){
x.spatialperm.empiricaly %>% group_by(as.character(p)) %>% #group by permuted S-A axis bins
do(mean.tvalue.significant = mean(.$empirical.y, na.rm = TRUE)) %>% #calculate mean t-value
unnest(cols = c("mean.tvalue.significant")) %>% t()}) %>%
t() %>% as.data.frame()
envSES_effects_SAaxisbins_permutedx <- envSES_effects_SAaxisbins_permutedx %>% mutate_if(is.character, as.numeric) #make it numeric
envSES_effects_SAaxisbins_permutedx <- envSES_effects_SAaxisbins_permutedx[ , !c(TRUE,FALSE) ]
colnames(envSES_effects_SAaxisbins_permutedx) <- c("SA.bin1.tvalue","SA.bin2.tvalue", "SA.bin3.tvalue", "SA.bin4.tvalue", "SA.bin5.tvalue")
## compute mean t-value for each S-A axis bin using spatially permuted y inputs
### set up spun (y) and empirical (x) df
y.spatialperm <- as.data.frame(y.perm) %>% set_names(sprintf("perm%s", seq(from = 1, to = ncol(y.perm))))
y.spatialperm.empiricalx <- cbind(x, y.spatialperm) %>% as.data.frame()
colnames(y.spatialperm.empiricalx)[1] <- c("empirical.x")
### compute mean t-values for S-A axis bins for all spins
envSES_effects_SAaxisbins_permutedy <- y.spatialperm.empiricalx %>% group_by(empirical.x) %>% #group by empirical S-A axis bins
dplyr::summarize(across(everything(), \(x) mean(x, na.rm = TRUE))) %>% #calculate mean t-value for permuted data
select(-empirical.x) %>% t() %>% as.data.frame() %>%
set_names(c("SA.bin1.tvalue","SA.bin2.tvalue", "SA.bin3.tvalue", "SA.bin4.tvalue", "SA.bin5.tvalue"))
histogram.data <- rbind(envSES_effects_SAaxisbins_permutedx, envSES_effects_SAaxisbins_permutedy)
saveRDS(histogram.data, "/cbica/projects/thalamocortical_development/thalamocortical_results/environment_results/envtvalue.SAbins.spatialnulls.HBN.rds")histogram.data <- readRDS("/cbica/projects/thalamocortical_development/thalamocortical_results/environment_results/envtvalue.SAbins.spatialnulls.HBN.rds")
envSES_effects_SAaxisbins_permutedx <- histogram.data[1:10000,]
envSES_effects_SAaxisbins_permutedy <- histogram.data[10001:20000,]
nperm = dim(perm.id.full)[2]#Function to plot the null distribution of environment effect t-values for a given S-A axis bin and test whether the true t-value is significantly smaller or greater in magnitude than the null (permutation based p)
enveffects_SAbin_enrichment <- function(quintile, enrichment_type){
SAbin.tvalue.empirical <- envSES_effects_SAaxisbins %>% filter(SA.axis.bin == quintile) %>% select(mean.tvalue.significant)
#Histogram
SA.colorvector <- c("#FFC228", "#F9D58C", "#e9dcf2", "#BE93C4", "#6F1382")
null.tvalue.histogram <- histogram.data %>% select(sprintf("SA.bin%s.tvalue", quintile)) %>%
ggplot(aes(x = .[,1])) + geom_histogram(fill = c(SA.colorvector)[quintile]) +
theme_classic() +
geom_vline(xintercept = SAbin.tvalue.empirical$mean.tvalue.significant, linewidth = 0.25) +
xlab("Environment effect (t-value)") +
ylab("Number of nulls") +
scale_x_continuous(breaks = c(2, 3, 4), limits = c(2, 4)) +
theme(
legend.position = "none",
axis.text.x = element_text(size = 5, family = "Arial", color = c("black")),
axis.text.y = element_blank(),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line.x = element_line(linewidth = .2),
axis.line.y = element_blank(),
axis.ticks = element_line(linewidth = .2))
#Permutation-based p
if(enrichment_type == "smaller"){
p.perm.xy <- sum(envSES_effects_SAaxisbins_permutedy[,quintile] < SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
p.perm.yx <- sum(envSES_effects_SAaxisbins_permutedx[,quintile] < SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
perm.p <- (p.perm.xy+p.perm.yx)/2
}
if(enrichment_type == "larger"){
p.perm.xy <- sum(envSES_effects_SAaxisbins_permutedy[,quintile] > SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
perm.p <- p.perm.xy
p.perm.yx <- sum(envSES_effects_SAaxisbins_permutedx[,quintile] > SAbin.tvalue.empirical$mean.tvalue.significant)/nperm
perm.p <- (p.perm.xy+p.perm.yx)/2
}
permutation.results <- list(null.tvalue.histogram, perm.p)
return(permutation.results)
}S-A quintile 1
## [[1]]
##
## [[2]]
## [1] 0.19655
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/quintile1_enveffect_spins_hbn.pdf", device = "pdf", dpi = 500, width = 1.1, height =.8)S-A quintile 2
## [[1]]
##
## [[2]]
## [1] 0.1787
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/quintile2_enveffect_spins_hbn.pdf", device = "pdf", dpi = 500, width = 1.1, height =.8)S-A quintile 3
## [[1]]
##
## [[2]]
## [1] 0.19985
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/quintile3_enveffect_spins_hbn.pdf", device = "pdf", dpi = 500, width = 1.1, height =.8)S-A quintile 4
## [[1]]
##
## [[2]]
## [1] 0.5911
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure8/quintile4_enveffect_spins_hbn.pdf", device = "pdf", dpi = 500, width = 1.1, height =.8)S-A quintile 5
## [[1]]
##
## [[2]]
## [1] 0.03145